EPJ manuscript No. 

(will be inserted by the editor) 



The one-dimensional Lennard-Jones system: 

collective fluctuations and breakdown of hydrodynamics 



Stefano Lepri a , Paolo Sandri, and Antonio Politi 

Istituto dei Sistemi Complessi, Consiglio Nazionale delle Ricerche, Sez. Territorial di Firenze 
and Istituto Nazionale di Ottica Applicata, largo E. Fermi 6 1-50125 Firenze, Italy 

Received: February 2, 2008 

Abstract. The dynamical correlations of a model consisting of particles constrained on the line and interact- 
ing with a nearest-neighbour Lennard-Jones potential are computed by molecular-dynamics simulations. 
A drastic qualitative change of the spectral shape, from a phonon-like to a diffusive form, is observed upon 
reducing the particle density even ad moderate temperatures. The latter scenario is due to the sponta- 
neus fragmentation of the crystal-like structure into an ensemble of "clusters" colliding among themselves. 
In both cases, the spectral linewidths do not follow the usual q 2 behaviour for small wavenumbers q, 
thus signalling a breakdown of linearized hydrodynamics. This anomaly is traced back by the presence of 
correlations due to the reduced dimensionality. 

PACS. 05.60.-k Transport processes - 66.10.Cb Diffusion and thermal diffusion 

1 Introduction to vanish for a finite one (see and references therein). 

Another example is the divergence of the thermal con- 
Relaxation and transport phenomena in reduced spatial ductivity coefficient with the length observed in chains 
dimension (D < 3) are often qualitatively different from of anharmonic oscillators an d hard-points gas [TU], 
their three-dimensional counterparts. For concreteness, imag-The same type of anomaly has been detected also for a 
ine a large set of impenetrable spheres confined within a quasi-l£> model consisting of spheres confined in a nar- 
narrow channel. If the mutual passage of particles is for- row channel 

bidden, the motion of the spheres is necessarily correlated, Computer simulations of simple toy models is an in- 
even at long times, because the displacement of a given va i uab i e way of attacking those problems. In particular, 
particle over a long distance necessitates the motion of one would like to understand the conditions under which 
many other particles in the same direction. This is a doc- those anoma iies occur and possibly to classify the possi- 
umented effect, for example, in single-filing systems where ble un i ver sal features. In this paper we consider a simple 
particle diffusion does not follow Fick's law [HE]. Another model of point part icles constrained on the line and inter- 
related instance is the enhancement of vibrational energy acting with their nearest neighbours through a Lennard- 
transmission in quasi-lD systems like polymers or in- Jones force This type of phenomenological interaction has 
dividual carbon nanotubes |3j . been throughly studied for decades by molecular-dynamics 
The distinguished signature of those effects is in the me thods \T2\. However its one-dimensional version has re- 
long-time behavior of the associated correlation functions ceived little attention so far. Previous studies focused on 
0. As it is known, the latter may display long-time tails an harmonic effects [HHH| and transport properties . 
leading to ill-defined transport coefficients or, more gen- In this reS pect, some evidence that the energy current au- 
erally, to the breakdown of customary hydrodynamics. In- tocorrelation (the Green-Kubo integrand) shows a long- 
deed, power-law decay of correlations is expected to be a time tail of the above mentioned type has been provided 
generic feature of one-dimensional systems in presence of t oo 16 . The same system has been also proposed as a toy 
conservation laws 0. One important consequence of such model to desC ribe fracture nucleation [rflfT^j . 

long-ranged correlations is that physical properties may _„ . , . , „ „ . , , . „ 

■ -c j j j-u 4. • i ,i . .-I Of particular relevance m what follows is the work of 

significantly depend on the system size and that the tner- _. , 1 . ... , . . , . ... 

j . j • n •■ v i tJishop and collaborators 19,20. At variance with our 
modynamic and mhmte-time limits may not commute. , , , . , , , — ^~ . , . , , 
t— i • ,, j ,. , j.rr . rc ■ ■ model, they considered the case m which the interaction 
tor instance, the tagged-particle diffusion coefficient m ' f „ . , . _ , „ l . ... 
n , . a ., r . n • r j- j is extended to all particle pairs. 1 hey noticed that the litc- 
u = 1, that is finite tor the infinite system is found instead ; . r , , t, . J . . .,, 
; time of long-waveiength fluctuations does not scale with 

a Electronic address: stefano.lepri@isc.cnr.it their wavenumber q as q~ 2 , as expected from standard 
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hydrodynamics, but rather as a nontrivial power g _/i . We 
anticipate that the results presented henceforth are qual- 
itatively consistent with theirs. However, the estimate of 
the exponent fi = 1/3 given in Ref. (20] is significantly 
different from the value found here, fi ~ 1.5, which is con- 
sistent with previous measurements in chains of coupled 
oscillators EHE21 This issue is relevant especially for as- 
sessments about the universality of the scaling laws. 

The present paper is organized as follows. In section 2 
we define the model and its physical parameters. Section 3 
is devoted to a discussion of the nature of the ground state 
at zero temparature. This is important to understand the 
effect of density changes on the system dynamics. In Sec- 
tion 4 we present our simulation results and show how 
the low-density properties are related to the kinetics of 
particle clusters. Finally, we summarize our results in the 
concluding Section. 

2 The model 

We consider an array of N point-like identical atoms or- 
dered along a line. The position of the n-th atom is de- 
noted with x n . By fixing the mass, without loss of gen- 
erality, equal to unity and assuming that interactions are 
restricted to nearest-neighbour pairs, the equations of mo- 
tion write 

x n = -F n + F n -i ; F n = -V' (x n+1 - x n ) , (1) 

where V'(z) is a shorthand notation for the first derivative 
of the the interparticle potential V with respect to z. The 
particles are confined in a simulation "box" of length L 
with periodic boundary conditions 

x n +N = x n + L . (2) 

Accordingly, the particle density d = N / L is a state vari- 
able to be considered together with the specific energy 
(energy per particle) that will be denoted by e. 

In the present work we focus on the Lennard-Jones 
potential that in our units reads 

For computational purposes, the coupling parameters have 
been fixed in such a way as to yield the simplest form for 
the force. With this choice , V has a minimum in z = 1 
and the resulting dissociation energy is Vq = 1/12. No- 
tice that for convenience we set the zero of the potential 
energy in z — 1. The presence of the repulsive term in 
one dimension ensures that the ordering is preserved (the 
particles do not cross each other). 

Before closing this section, let us recall that most pre- 
vious studies |16U15U19II20| dealt with the case in which 
the interaction is not restricted to nearest neighbours but 
rather extends all particle pairs. In practice, as pointed 
out in Refs. |19U2(J| . even in this case the interaction is 
limited to about 2-3 neighbours. Therefore, we do not ex- 
pect that the results reported below will be significantly 
altered when taking into account the interaction among 
all pairs. 



3 The ground state at T = 



In order to understand the physical features of the model 
let us briefly discuss its equilibrium properties. Since the 
system is one dimensional with only nearest-neighbour 
and short-ranged interaction, no actual phase transition at 
finite temperature, T > 0, can occur. Nevertheless, we will 
see that the dynamics of collective modes may strongly 
depend on the state variables. Therefore, even if it is not 
legitimate to speak of "phases" in a strict thermodynamic 
sense, it is sensible to distinguish between regions where 
the dynamics is qualitatively different. 

Let us start discussing the equilibrium configurations 
at T = as a function of the particle density d. This issue, 
which is important for the choice of the initial conditions 
in microcanonical simulations (see below), has been dis- 
cussed in Ref. for the Lennard-Jones chain where the 
pair interaction is among all particles. In view of the short- 
range nature of the potential we do not find significative 
difference with the case at hand here. For convenience, 
we summarize the main results in Fig^ Three density 
regimes are distinguished: 

— For d > 1 the ground state consists of equally spaced 
particles at a relative distance a = l/d (homogeneous 
solution). The chain is under compression and the total 
energy decreases upon decreasing d up to the minimum 
value which is attained for d = 1. 

— For < d < 1 the homogeneous solution (chain un- 
der tension) becomes a relative minimum. The ground 
state consists of equally spaced particles ad a distance 
approximatively a rs 1 except for a couple which lies 
at a large distance. In other words, the minimal en- 
ergy configuration is attained by breaking the chain at 

one bond. At the critical value d = d* = \J~^g — 

0.901971 . . ., the homogeneous configuration undergoes 
an instability and disappears. 

— For d < d* the broken chain state is the unique mini- 
mal energy configuration. 

Notice that for d < 1 the ground-state energy is basically 
up to terms 0(1/N) since only a couple of particles out 
of N sits at a distance which is significantly larger than 
1. Furthermore, as noted in Ref. |23j . further fragmenta- 
tion in two or more shorter chains cannot produce further 
energy minima. 

The phase diagram of Fig. suggests a first-order 
phase transition at d = 1 where the second derivative 
of the energy density with respect to d undergoes a jump 
discontinuity. As this quantity should be proportional to 
the Young's modulus this physically means that the chain 
undergoes a fracture and looses its elastic tension. 

At finite temperature we expect that these features 
should be washed out but, as observed for other mod- 
els [21], the remnants of the transition should somehow 
manifest themselves in the dynamics. 
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4.1 The d = 1 case 
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Fig. 1. The T = minimal-energy configurations of the 
Lennard-Jones chain as a function of the particle density. 



4 Dynamical correlations 



We have performed equilibrium microcanonical simula- 
tions by integrating Eqs. Q (with periodic boundary con- 
ditions (J2J ) by means of a fourth-order symplectic al- 
gorithm For the different values of density consid- 
ered henceforth, initial conditions were chosen to be in 
the ground state described in the previous section. The 
inital velocities were drawn at random from a Gaussian 
distribution and rescaled by suitable factors to assign the 
kinetic energy to the desired value and to set the total 
initial momentum equal to zero. A suitable transient is 
elapsed before acquisition of statistical averages. Conser- 
vation of energy and momentum was monitored during 
each run. This check is particularly crucial at high en- 
ergies/densities where the strongly repulsive part of the 
force comes into play and may lead to significant inaccu- 
racies. The chosen time-step (0.005-0.05) ensures energy 
conservation up to a few parts per million in the worst 
case. 

We computed the dynamical structure factor, namely 
the square modulus of temporal Fourier transform of the 
particle density 



which is defined as 



77 expH^n) 



S(q,oj) = (\p(q,uj)\ 2 ) 



(4) 



(5) 



The square brackets denote an average over a set of inde- 
pendent molecular-dynamics runs. By virtue of the peri- 
odic boundaries, the allowed values of the wavenumber q 
are integer multiples of 2n/L. The reliability of the spec- 
tra have been checked against different choices of the run 
durations and sampling times. In some particular cases 
we also verified that the results are not affected by data 
windowing that, in principle, may affect the measured 
linewidths. 



In this case we expect a "lattice-like" behaviour with each 
particle oscillating around its equilibrium positions. In 
other words, we can introduce the change of coordinates 
x n = u n + na where a = 1/d is the constant lattice spac- 
ing. The computation of S(q,u) can be more conveniently 
performed by resorting to the collective coordinates 



U(q) 



1 N 

— = ^2 u n exp(-iqn) q = 

V n—1 



2vrfc 
~Na 



(6) 



for k being an integer comprised between —N/2 + 1 and 
N/2. This latter expression is computationally more effi- 
cient than the definition Q. In fact, standard Fast Fourier 
Transform routine can be used to evaluate U(q) (provided 
that N is a power of 2). On the other hand, by expand- 
ing to the leading order in q, it is seen that the structure 
function (J5J is proportional to q 2 {\U{q)\ 2 ). We checked nu- 
merically that this approximation is very accurate in the 
range of q values considered henceforth. 

In Fig. we report two representative data sets for 
low (e = 0.02) and high energies (e = 0.2) compared with 
the well depth. The most distinguished feature of the the 
spectra is a narrow phonon peak. For the low energy case, 
its frequency fi agrees with the one computed by the har- 
monic approximation of the Lennard-Jones potential: 



Q{q) = 2y/V"(l) I sin 



qa , 



(7) 



Indeed, the estimated sound speed c = 2.77 ± 0.06 is 
only slightly larger the value obtained from this latter 
expression (c = y/6 w 2.49...). This indicated that non- 
linear terms are weak in this energy range. At higher en- 
ergies, the particles are less sensitive to the bottom part of 
the potential well and the phonon peaks shifts to higher 
frequency. The data in Fig. |2h correspond to e = 0.2, 
more than 20 times the well depth; comparing the values 
of Fig. it can be ascertained that the frequencies are 
roughly 50 % larger than in the previous case. 

For small enough wavenumbers, besides the phonon 
peak, a small zero-frequency component appears whose 
relative weight increases at larger energies (see again Fig.EJ. 
This feature, that has not been reported in previous works 
on the Lennard-Jones chain [2| , is reminiscent of the well- 
known triplet structure observed in fluids where a cen- 
tral Rayleigh peak is accompanied by two narrow Bril- 
louin lines placed at ±cg [201 • Moreover, according to 
the standard hydrodynamic theory, the ratio of the ar- 
eas under half the Rayleigh peak and one Brillouin peak 
is C p /Cv — 1. On the other hand, the specific heats ratio 
is expected to be very close to unity for a "crystal-like" 
structure as the one we are facing in our simulation. This 
is thus consistent with the intuitive idea that a sizeable 
central component should occur only when the systems is 
"fluid enough" i.e. when fluctuations of the particles' po- 
sitions are large. This is expected for large temperatures 
and/or small densities . In the next section we will see 

The existence of a central component also means that there 
is a coupling between density and energy fluctuations which are 
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that this expectation is further confirmed by the simula- 
tions. 




Fig. 2. Structure factors for d = 1, N = 8192 and different qs 
corresponding to indices k = 1,2,4,8 in eq. @ (left to right). 
Microcanonical simulations are performed for the energy den- 
sities e = 0.02 (a) and e = 0.2 (b). The spectra are averaged 
over an ensemble of about 100 initial conditions. 



In Fig.|3we report the lincwidths of the Brillouin peaks 
as a function of the wavenumber q and for the two en- 
ergy values chosen above. The data reported there corre- 
spond to those wavenumbers for which the linewithds are 
small enough with respect to the peak frequency (typically 
less the one tenth) to allow for a meanigful estimate. The 
linewidths are computed by evaluating the frequencies at 
which the spectrum is one half of its maximal value. In the 
cases in which the spectral resolution is not high enough 
we used by a Lorenzian fit to improve the accuracy. Inci- 
dentally, we noticed that the the fitting is reasonably good 
only in the peak region while substantial deviations on the 
tails are observed. 

Remarkably, the linewidths do not scale as q 2 as ex- 
pected from standard hydrodynamics but rather as q^. 
For e = 0.02 our best-fit estimate is fi = 1.47 ± 0.05. To 



better appreciate the reliability of this value, in the inset 
of Fig. |3 we report the logarithmic derivative of the data 
that show a plateau around 1.5. The data for e = 0.2 are 
instead less conclusive. Although a power law fit still yields 
a comparable value (p — 1.41 ± 0.07) the plot displays a 
residual curvature that indicate some relevant sublcading 
corrections. 




the hydrodynamic fields of our model. To support this inter- 
pretation we computed also the spectra of the local energy field 
in some cases. They do display a similar structure with a large 
component at the phonon frequency. 



Fig. 3. The linewidths of the phonon peak for d = 1, e = 0.02 
and e = 0.2. Error bars are reported only when significantly 
larger than the symbols. The inset shows the logarithimc 
derivative evaluated by finite differences. 



4.2 The d < 1 case 

At finite temperature, when the density is lowered be- 
low d = 1, the periodic ground state destabilizes and 
this leads to a completely different dynamics. To avoid 
the additional features connected with the presence of 
the metastable branch (see Fig. ^) we consider the case 
d < d*. Fig. 01 shows the evolution of the particle density 
field starting from the minimal-energy configuration for 
d = 0.8. After a transient, the system spontaneously frag- 
ments in a series of clusters of different lengths. Within 
each cluster, the particles are separated by an average 
distance equal to the equilibrium distance of the Lennard- 
Jones potential (1 in our units). 

The relaxation is thus ruled by the clusters' kinetics. 
The latter consists of two different processes, fragmenta- 
tion and collisions whose characteristic times we denote 
by tf and tc respectively. Let us start giving an estimate 
of the average number m of particles in each cluster. For 
convenience, we write m = l/f where / is the average 
fraction of broken bonds. The average time between col- 
lisions is given by the ratio between the typical clusters' 
separation and the typical speed v. A straightforward cal- 
culation yields 

(remember that the equilibrium distance is 1 in our units). 
On the other hand, we expect fragmentation to be a ther- 
mally activated process with an activation energy equal 
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Fig. 4. The evolution of the particle density for d = 0.8, N = 
512 and e = 0.02. Along the horizontal direction a black pixel 
is drawn at each particle location. Time increases downwards. 



to the well depth Vq. Accordingly, 
t~f = rf exp (3V 



where r is a suitable prefactor setting the typical "attempt 
time" for the process and (3 is the inverse temperature. 
The condition for kinetic equilibrium is obtained by letting 
t~c = Tf, 



/ = 



1 - d 



exp 



2 



(10) 



Notice that, as a consequence of this estimate, the charac- 
teristic time is proportional to y/l/d — 1 exp(/3Vb/2) and 
that thermalization may become very slow for low tem- 
peratures and densities. 

This prediction is in fairly good agreement with the nu- 
merical data (Fig. EJ). There, we plot the average fraction 
of broken bonds / as a function of the inverse kinetic tem- 
perature. The quantity / is measured by counting at each 
time the number of pairs whose distance is larger than 
some prescribed threshold that we fixed equal to 1.5 (some 
50% above the inflection point of the Lennard- Jones po- 
tential). From the above discussion, it is clear that the 
choice of well-thermalized initial conditions is crucial. In- 
deed we found that the time to reach the equilibrium value 
of / increases (about linearly) with the number of parti- 
cles. In Fig.[S^ it is shown that / is an extensive parameter 
whose equilibrium value can be evaluated already for sys- 
tems of a few hundred particles. In Fig. 0J> we check that 
the scaling behaviour predicted by the kinetic argument 
reported above is in agreement with the simulation data. 
Notice in particular that the factor 1/2 in the exponential 
term is very well accounted for by the data and that the 
prefactors display a rather weak dependence on (3. 

In view of the above dynamical features, we expect 
that the the spectra of density fluctuations should be qual- 
itatively different from the d = 1 case. This is confirmed 
by the simulation reported in Figs. El (in this low-density 
regime we use the definition JSJ of S(q,uj)). Indeed, even 
at low energies, a large central component appears which 
we associate with the diffusive behaviour of the clusters. 
The secondary Brillouin peaks occur only for small enough 
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Fig. 5. The average fraction of broken bonds / as a function 
of the inverse kinetic temperature /3; (a) dependence of / on 
the number of particles N for fixed density d — 0.5; (b) check 
of the formula (TJ3 for fixed N = 256. 



wavenumbers and are weaker: for the smallest wavenum- 
ber shown in Fig. the secondary peak has about 40% 
spectral power of the central component. This observa- 
tion is consistent with the intuitive expectation that at 
lower densities the system should respond like a fluid. 

The behavior for e = 0.2 is shown in Fig.Et>. The struc- 
ture factor is qualitatively similar to the case of Fig. [2b- 
This is expected, since in the high energy limit the model 
approach the case of hard points colliding elastically. It is 
thus plausible that the dynamics is practically indepen- 
dendent on the particle density. Notice also that in this 
limit the model becomes an almost-integrable dynamical 
system and some some pathological behaviour (e.g. slow 
thermalization) may be expected. 

In analogy with the case d = 1 we analysed the de- 
pendence of the spectral linewidths on the wavenumber, 
considering only those q values that correspond to narrow 
enough lines. Fig0 shows the half-widths of the central 
(Rayleigh) peak as a function of the wavenumber q for 
the e — 0.02. In the accessible range we find again a non- 
trivial scaling q^ with ji = 1.60 ± 0.05. 
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Fig. 6. Structure factors for d = 0.8 and different qs cor- 
responding to indices k — 1,2,4,8 in eq. (left to right). 
Microcanonical simulations are performed for the energy den- 
sities e = 0.02, N = 1024 (a) and e = 0.2, N = 4096 (b). The 
spectra are averaged over an ensemble of about 500 and 100 
initial conditions respectively. 





Fig. 7. The linewidths of the central peak for d = 0.8 and 
e = 0.02. Error bars are reported only when significantly larger 
than the symbols. 



5 Conclusions 

Our numerical analysis of the one-dimensional Lennard- 
Jones system has shown that long-ranged correlations in- 
duced by the reduced spatial dimensionality may signif- 



icantly affect the density fluctuations. The most striking 
feature is the anomalous scaling of the Rayleigh and Bril- 
louin peak widths in the hydrodynamic limit, q — > 0. The 
physical meaning of this behaviour is understood by re- 
calling that, in the standard framework, the linewidths 
are connected to transport coefficients For instance 
the width of the Brillouin peak is Fq 2 where F is the 
sound attenuation constant. The anomalous scaling can 
be recasted in terms of a wavenumber-dependent constant 
r(q) ~ q^~ 2 . Since our simulations clearly indicate that 
fi < 2, this implies that r diverges in the q — > limit. 
In other words, a long- wavelength disturbance is damped 
on a typical distance which becomes very large. In some 
sense, one may think of this as a superdiffusive process, in- 
termediate between standard diffusive and ballistic prop- 
agation. Our result is thus closely related to the analysis 
performed in Ref. [27| for the hard-point gas. 

This feature has been previously observed also in other 
one-dimensional lattice models |21U22| . A similar anoma- 
lous behaviour of the viscosity of a ID lattice gas has 
also been reported in Ref. [2Hj- However, the exponent \x 
found in the present work is about 10% smaller then those 
measured in Refs. |21U22| (1.5 against 1.67). Remarkably, 
both values are very close to theoretical estimates pre- 
sented previously in the literature. Ernst claimed that 
/i = 5/3, while the more refined mode-coupling analysis 
by Wang and Li gave (x = 3/2 for the specific model 
they considered. We anticipate that this latter value is 
actually supported by numerical solution of the mode- 
coupling equations |31| that, according to Ref. |14j . should 
describe the dynamical properties of our model. This es- 
timate is in excellent agreement with the data reported 
above. Whether the small differences in the exponents 
are due to numerical errors and subleading corrections or 
they indicate the existence of two different "universality 
classes" is still an open question. 

For d — 1 the system retains the dynamical features 
of a one-dimensional crystal. Of course, by virtue of the 
short-ranged interactions, the model cannot have a gen- 
uine solid phase (Landau-Peierls instability). Nonetheless, 
we have found that for a finite system we can still reason 
in terms of an effective phonon dispersion and damping. 
In the low-density region {d < 1) we have shown how the 
relevant time scales are instead dictated by the collision 
and recombination processes of one-dimensional "clusters" 
of particles. Such processes may lead to a remarkable in- 
crease of the thermalization times when the kinetic tem- 
perature becomes smaller than the binding energy (see 
Eqs. EJandEJ). 

In the present study we limited ourselves to the case in 
which the system is initialized close to the ground state. 
It would be interesting to examine the relaxation from the 
metastable state which exists for e?* < d < 1. In this re- 
spect there is a connection with the work of Oliveira [T%] 
that studied the fracture nucleation in the same model. At 
variance with what discussed here, he starts from the uni- 
formly stretched chain, namely from the metastable state 
illustrated in the phase diagram of Fig.Q] He founds that 
the breaking time is orders of magnitude longer than what 
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expected from Kramers-type estimate. This is presumably 
to be traced back to subtle long-ranged correlations of the 
very same type studied here. 

To conclude, we wish to mention the fact that the 
anomalous scaling may manifest also at the level of en- 
ergy transport in this model. Indeed, nonequilibrium sim- 
ulations show that the thermal conductivity coefficient 
diverges with the system length also in the low-density 
regime 32 . This is a further evidence that nonequilib- 
rium processes in one dimension are peculiar and deserve 
a special attention. 
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